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INTRODUCTION 


The common practice in applying pseudospectral methods to partial dif- 
ferential equations is to satisfy the equation at the interior nodes and to 
impose the boundary conditions at the boundary* This procedure does not take 
into consideration that the differential equation is satisfied at points arbi- 
trarily close to the boundary. In [4], the authors discussed the advantages 
of imposing a combination of boundary conditions and the equation itself at 
the boundary nodes, for Chebyshev approximations of the Laplace equation with 
Neumann conditions. Here we analyze the same idea applied to the linear 
hyperbolic equation 

! u t = u x l x l £ t > o, 

u(x,0) = f(x), 
u(l,t) = g(t). 

We assume that the collocation points are the Gauss-Lobatto Chebyshev 
quadrature nodes. The stability of the method, with the commonly used 
boundary treatment, i.e., imposing u^(l,t) = g(t) was analyzed in [10]. 
Here we show the convergence of the method for the new boundary treatment, 
namely: 

TT (1 » t) " JT (1 ’ c) + ^V 1 * 0 “ g(t)) = 0 

where o is positive and large enough. A preliminary theoretical discus- 
sion in Section I and numerical experiments in the last Section show the 
effectiveness of the method. In Section III we use the results obtained for 
the hyperbolic equations to show that the second derivative matrices, cor- 
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responding to the Neumann conditions with the new approach, have real and 
negative eigenvalues. The analogous result for the classical way to impose 
boundary conditions was previously proven in [8]. 


1. DESCRIPTION OF THE NEW METHOD 

In order to illustrate the new method of imposing boundary conditions and 
to explain what gain can be realized by this method, we begin with the 
following simple time independent problem 


( 1 . 1 ) 


j°x ■ f 

(u(l) = 0, 


where f€C s ([-l,l]) is given (s _> 0). 

In the standard pseudospectral Chebyshev method (see for instance [6]), 
we look for a polynomial of degree N, say vjq, such that 

(Xj) j = !»••• ,N, 




dv N r i 
IT ( V ' 


((a) 

(1.2) 

1(b) 

y» - 0. 

where x. 

3 

= C0S N^ ’ 

j = 0,1, ••• 

in [-1,1]. 

In order to 

determine 


unknown point values v N (xj) using the Lagrange interpolation polynomial 


N 


V N <X> ‘ S V’ScV’ 0 ' 

k=0 


where 


g k ( x) = 


-(-l) k (l - x 2 )T'(x) 
c fc N 2 (x “ x k ) 
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t 


[ 


Here is the N-degree Chebyshev polynomial and Cj = 1 if 1 j £ N-l, 

while c 0 = c N = 2. Therefore 


dV N ? d % 


2 W *r ( V 

c=0 J 


j = 1,« 


,N. 


Upon substituting the above relations in (1.2), we get a linear system of 
equations for the point values Vjj(x^). We note that in (1.1) the dif- 
ferential equation holds in any arbitrary neighborhood of the boundary, where- 
as in (1.2) the differential equation is satisfied not further than x = Xj. 
We did not require, for instance, that the equation could also be satisfied 

at xq = 1. We propose now another procedure that takes into account the 
*- 

differential equation at the boundary as well as the boundary condition# 

In our new method, we seek an N-degree polynomial u^ such that 


i 

I (1.3) 


(a) 

du » , , 

= £(x ) 

j = !,••• ,N 

(b) 

du N 

sr (I) ■ 

- «v u = 

f(l), 


where a > 0 is a suitable constant depending on N, to be determined 

1 du N 

later. By writing the equality (1.3)(b) as — (— f ) ( 1 ) = u N (l), we note 

that (1.2) is obtained from (1.3) by letting a +°>>. We remark that the 
solution of (1.3) satisfies neither the boundary condition nor the equation 
at x » 1; if the method converges both will be satisfied as N + +». 

To show the advantage of the new procedure, we give in Figure 1.1 the 


plot of the error 
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(1.4) E » E(a) = (|r I (U - u ) 2 (x ) -^) 1/2 

.1=1 3 j 

multiplied by 10 5 versus a, for f(x) = sin(x - 1) and N = 8. The 

point Xq is not taken into consideration in the sum because the exact 

solution is known there. It is clear from the figure that E(a) is not 

monotone in a and there exists a = a min which minimizes E. In 

particular we have E ( a m i n ) < E(+“). Further experiments indicate that, in 

o 

terms of N, ct . increases like N . 

min 



Figure 1.1 - Behavior of the error versus a. 

We would like to explain why the procedure (1.3) should be, in general, better 
than (1.2). We start by noting that if f is a polynomial of degree N-l at 
most, then both (1.2)(a) and (1.3)(a) hold, not only at the grid points x ^ , 


but for every x since both sides of the equations are polynomials of de- 

^ U N 

gree N-l. In particular (1) = f(l), thus by (1.3)(b) we get u N (l) “ 0» 

is a 


up to 
poly- 
solu- 


It is clear that the solution of (1.2) is 

(1.7) v jj(x) = 0, ¥x. 

On the other hand, the solution of (1.3)(a) is a constant and from (1.3)(b) we 
get 

(1.8) u N (x) = - | , ¥x. 

Set — = 8 , then the error is given now by 
a ’ 

i ■ ( j I* <«*,> + e> 2 r-) 1/2 - 
B j-i 3 j 


leading to the conclusion that u^(x) = v^(x), ¥x. Suppose now that f 
polynomial of degree N. We can assume, because of the linearity, that 


(1.5) 


f(x) 


(1 + x)T'(x) 
2 

2ir 


Hence: f (x_. ) = 0, j ■ 1,««»,N and f ( 1 ) = 1. Any other polynomial, 

a constant factor, can be obtained from (1.5) by adding some suitable 
nomial of lower degree. In this case it is easily verified that the 
tion U of (1.1) is given by 


(1.6) U(x) = 2 [ 2 N+1 ^ N+1 


2N' 


T„. , (x) + T n (x) +1^ T n-1 


, \ 2N 2 -1-| 

(x) - -g j . 

N-l 


(1.9) 
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To minimize E one has to choose 3 as the negative mean of U, namely 


min 


?r £ U(x.) — 

N i-i 3 c i 

N 

ir r 

N I c -j 
j*l 3 


N 


I U(x)f 
_ j-o 3 i 


N 


I f 

3=1 j 


and an easy calculation shows that 


(1.10) 


min 3 


min 


^ (N 2 -l) 
2N -1 


N 


This explains the behavior of a m i n as a function of N. 

In Table 1.1, we summarize the results of another experiment. This time 
3 1/2 

we took f(x) = - (1 “ x) , with the boundary condition U(l) = 1, so 

that the solution was U(x) - (1 - x )^/2 + i. We have tried the two different 
ways of imposing boundary conditions, i.e. 


(1.11) v N (l) = 1 


d u N 

(1*12) dx “min^N^ 1 ^ “ O = f(l), (where o min is given by (1.10)), 


and we varied the number of grid points N 


N 

Condition (1.11) 

Condition (1.12) 

2 

0.281837 

0.243315 

4 

0.338991E-01 

0. 183894E-01 

6 

0. 995738E-02 

0.478991E-02 

8 

0.418410E-02 

0. 190129E-02 


Table 1.1 - Comparison of the errors between the two ways of imposing boundary 
conditions. 


2. THE TIME DEPENDENT PROBLEM 

In this section, we show how to apply the new procedure of setting the 
boundary conditions, described in the previous section, to a hyperbolic equa- 
tion. An analysis of the convergence of this method will be carried out for 
Chebyshev approximations. Consider the equation 

(2.1) ( U t = U x M < t > °> 

< U(l,t) = g(t), 

( U(x,0) = f (x) . 


The pseudospectral semi -discrete approximation to (2.1) we suggest in this 
paper involves seeking a polynomial Ujj of degree at most N such that 
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(2.4) = cos II- J - 0,1,... ,N, 

determine a different version (see [7]). Pseudospectral Legendre is defined 
by choosing Xj to be the extrema of the N-th degree Legendre polynomial. 

We would like to show here the convergence of the solution Ujj(x,t) of 
(2.2) to U(x,t) defined in (2.1) when N -► +°°, in the case of Chebyshev 
methods. The proof of the stability in the case a = +°° (i.e. , the common 

way of imposing boundary conditions) was discussed in [10]. Here we follow 
the same basic ideas. We give the proof with some detail, since it will be 
useful to further results concerning systems of differential equations. For 
the convenience of the reader, we also report the proof of convergence which 
never appeared in literature. We start with the following preliminary re- 


sults 
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Lenma 2.1: Let u N (x,t) be the solution of (2.2) when Xj are given by 

(2,3), then 


(2.5) 


3u N 3u N . (l+x)T'(x) 

TT = JT + T ^2 


where x = - a(u N (l,t) - g(t)). 


Proof : It is sufficient to note that (2.5) exactly coincides with (2.2), when 

evaluated at the collocation nodes. ■ 


We define now P^U as the polynomial interpolating U at the points 
cos , j = 0,1, Note in particular that (P^U)(l,t) = U(l,t), for 
any M. We are ready to write the error equation. 


Lemma 2.2: 


( 2 . 6 ) 


Let 


e N (x,t) 




u N (x,t) - P N _ 3 U(x,t), then 
(l+x)T'(x) 

■ o 5 e N (l,t) + Q(x, t) 

2N 


e N (x,0) 


P N f ~ P N-3 f ’ 


where Q(x,t) is a polynomial of degree N-3 ^iri x, given by 
= 9x ^ P N-3 U ^ “ P N-3 ( 97 U) * 


Proof : We apply %_3 to (2.1) to get 

a ( p n _ 3 u) a ( p n _ 3 u) 


(l+x)T'(x) 


2N 


(2.7) 


at 


ax 


- Q(x, t) - a[(P N _ 3 U)(l,t) - g(t)] 


y 
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with the initial condition t-^N- 3^^t-0 = ^N-3^* ^ ac ^-» no *- e that, since 

(P N _ 3 U)(l,t) = g(t), the last term that was introduced in (2.7) is zero. 
Hence, (2.6) follows from (2.7) and (2.5). ■ 

Next we will show that e N (x,t) tends to zero as N increases. The 
proof will be based on a careful energy estimate for (2.6). For this, we need 
the following lemmas. 

4N-1 

Lemma 2.3: Let w(x) = J b T (x), then 

k=0 


( 2 . 8 ) 


N 



“ ( V 

c . 

3 


f 1 w(x) 

- 1 fiT? 


dx + irb 2N , 


where c o = C N = ^ and - 1 f° r 0 < k < N. 

Proof: We test (2.8) for T^, k = 0,«**,4N— 1. If 0 _< k _< 2N— 1, (2.8) is 

a well-known quadrature formula (see [9]); if k = 2N, it is a trivial result 

by noticing that T 2 N (xj) = 1; if 2N+1 <k< 4N-1, then by writing T 2N+m = 

= 2T„T -Tom . (2.8) follows easily from the orthogonality of the Chebyshev 
N m ZN-m’ 

polynomials. I 

N 

Lemma 2.4: Let v(x) = £ a.T, (x), then 

k=0 

l^ 0 ^ (1 + V° - - 


(2.9) 
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(1 + x) ( 1 - $x)\ 


f 1 *- 

1 ~ X 


-dx +|[(1 


B ^ Na N 3 2 a N a N-l^’ 


for any 3 real. 


Proof: The result is an application of the previous lemma (see also [10]). 


Lemma 2.5: Let e^(x,t) be defined by (2.6). Suppose that 

— £ aiT-fx)) then 

k=0 


(2.10) 


dt ^ 2a N “ a N-l ^ ~ “4N(2a N - a^^), 


Pr °o f ; w e can argue as in [10] using (2.6) and the fact that Q(x,t) is a 
polynomial of degree N-3. | 


Lemma 2.6: Let e N (x,t) be defined by (2,6). then 


(2-11> ^ j=o V 1 + * 3)(1 ’ Sx J ) * H<x J ,t) ^ ‘V 0 ■ 

= (- 1 (l+x)(l-3x) 3e N a 2N-1 d (r) .2 . 

i e N 3x ^ i am At- (2 a v - a M _, ) + 


/l - 


T5n~ dt ' N “ a N-T 


- | n< 3B - 1 - |)aj. 


Proof: Combine the results of Lemma 2.4 and Lemma 2.5. 
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Theorem 2.1: Define 


.2 TT r 1 


(2.12) He N » N ^ d + X j^ x “ 2 x j^ e N^ x j ,t: ^ + " 16N ^ 2a N a N-P ’ 


and let 


K 1- z +z 
r mm 


_1 v ] 

Y K 2K 2., 1- 


where z m are the zeroes of and K Is chosen such that K _> N+l. 

Then we have 


(2.n) 2 dt «e N ! + 2 Y K^ e N ^ 1,<: ^ - 


<2 / (1 + x)(l - ^) 2 Q 2 (x,t)/l - x 2 dx, Vt > 0. 

-1 


Proof : We evaluate the equation in (2.6) at the points x ^ , then we multiply 

by (1 +Xj)(l - y Xj)e N (xj,t) and sum up over j = 0,.».,N to get 

( 2 .i4) - ^ c7 (1 + X j )(1 ” 2 x j )e N (x j ,t:) 3t e N (x j ,t:) 

= N E 77 ^ + X j^ x “ 2 X j^ e N^ X j ,<: ^ ITx e N^ X j ,t: ^ + 

J=0 j 

" 2N e N (1,t) + ¥ cT (1 + X j )(1 ” 2 x j )e N (Xj ,t:> . 


The right hand side of (2.14) is composed by three terms. We start by 

estimating the last term. First, we realize that the polynomial 
(1 + x) ( 1 - | x)e N Q is of degree 2N-1 and therefore, by Lemma 2.1, we have 
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F .l 0 cT (1 + x j )(1 ■ 1 x J )e N (x j ,t)Q(x j ,t) 


- / (1 + x)(l - 4*)e M (x,t)Q(x,t) — — . 

/TT? 

Upon using the Gauss quadrature formula based on z^, m = 1,...,K one gets 


l£ 


(l+x)(l - ■i-x)e„Q — — 

- 1 fTT7 


K . 

T l ( 1 + 0 ( 1 “ y z„)e m (z ,t)Q(z ,t) < 

*• m 2 m N m m — 

m-i « 


K 1-z +z . 
m m 2 , 


- 2K Ji 2(l-z ) e N (z m ,t) + W \ 7 “T (1+Z m )2(1 

m-I m m=l 1-z +z 


T z rt) Q (z_»t) < 
z m m — 


m m 


i5r X I, «--.><»v 2 « - 

m- X m m= 1 


k>V<v e > 


K 1-z +z 0 
mm2 


2K ^ 2n Z rT e N (z m >t) + 2 / ^ 1+x ^ X " y x ) 2 Q 2 (x,t) /l - 

m=l m - 1 


2 , 
x dx. 


For the first term in the right hand side of (2.14) we use the result in Lemma 
2.6 with 0 = X . Therefore, by ( 2 . 12 ) and the previous estimate, we get 


7 Ft le N ,2 -T / (l+x)(l - X x ) [e^(x,t) - e^(l,t)] 


dx 

/l-n 



1 v 2 ira 2 
2N ;a N “ "2F e N 


(l,t) 


+ 


(2.15) 
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+ w <1+x)<1 

m-1 m J* 


1 v2 2. * dx 

yx) Q (x,t) 


Integration by parts for the first term in the right hand side of (2.15) 
yields 


/ (l+x)(l - 


7° h. k N ( *' t) ■ s 0 ’ 01 t= 


1 f * 1-X+X 2 r 2, v 2 ( . v, dx 
= - 7/ [e N (x>t) " e N (1 » t)] — 


“ 2 ^ 2( 1-z ) e N (z m ,t:) Y K 1te N (1,t ^* 

m=l m 


where we noted that the last integrand is a polynomial of degree 2N+1 2K-1 

and therefore the Gauss quadrature formula is exact. Going back to (2.15), 
one finally gets (2.13). ■ 


Remark 2.1: It can be shown that H • II defined in (2.12) is actually a 


norm. In fact, it is possible to find a constant c, independent of N, such 


>_ c j e N (1 + x) 


'1 - x 


for every polynomial of degree at most N. 


Finally, by integrating (2.13), we get the main result of this section. 


Theorem 2.2: Let_ a be such that n(^ - Y^) > c » where c does not 

depend on N, then we have 
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(2.16) II e N ( • ,t)l 2 + C* f* e 2 (l,T)dx < 

t 1 ^ 

< IIP f - P -f| 2 + 4/ J (l+x)(l - y) 2 Q 2 (x,X )/l-X dxdx. ■ 

— N n-j 0-1 

The previous theorem is a convergence result by noting that the right hand 
side of (2.16) goes to zero in a spectral way (see for instance [1]). 

y n+i i 

Remark 2.2 : One can check that — converges to y when N goes to •+«. 

This means that, by taking a proportional to N^, the hypothesis of Theorem 
2.2 is satisfied. This assumption is similar to that made for the time 
independent problem (see (1.10)). H 


3. BOUNDARY CONDITIONS FOR ELLIPTIC EQUATIONS 

A theoretical analysis of the convergence for pseudospectral approxima- 
tions of the solution of Neumann problems, with a modified approach to treat 
the boundary conditions similar to that examined in the previous sections, has 
been developed in [4], Here we shall prove that the matrices relative to such 
approximations have real and strictly negative eigenvalues (note that, in the 
Chebyshev case, these matrices are not symmetric). To this purpose, we con- 
sider the parabolic equation 

(3.1) U t = U xx , |x| < 1 


with the Neumann boundary conditions 
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(3.2) 


U x (±l) = 0. 


The solution is determined up to a constant. The Chebyshev method with the 
new boundary treatment Involves seeking an N ttl degree polynomial u N such 
that 


(3.3) 


3u, 

at 


N 


a 2 u 


N 


ax 


at x = cos j = 1,... ,N-1, 


and 


(3.4) 


at 

at 


a x 2 


+ a = 0 


3x 

3u n 


* 2u n ” - 0 


at x = 1 


at x = -1, 


where a is a positive constant to be determined later on. The eigenvalue 
problems associated with (3.3) - (3.4) consist of finding a non-vanishing 
polynomial v, of degree at most N, such that 


(3.5) 


Xv = V xx at x = x j» J = 


and 


Xv - v + av = 0 


xx 


at x = 1 


(3.6) 


Xv - v - av = 0 

XX X 


at x = -1. 


The problem (3.5) admits the trivial solution X = 0. We will show that the 
other eigenvalues are real and strictly negative. We begin by noticing that 
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one can explicitly derive the characteristic polynomial of (3.5) - (3.6). In 
fact, (3.5) can be written as follows 

(3.7) Xv = v xx + aR + bS, a,b€lR, 


where 


R(x) = 


xT'(x) 


S(x) = 


t '(k> 

N 2 


Therefore, following [8], we have the next result. 


Lemma 3.1: The solution v of (3.7) is given by 


(3.8) v(x) = a p(x,y) + b q(x,y). 


where y 


1 

X 


and 



j 

p(x,y) 

r _(2k)> >. k+1 

= l R (x) v 


] 


k=0 

(3.9) 

j 

f 

q(x,y) 

- 1 S (2k> (x)/ +1 . 



V 

k=0 

Proof : 

We first note that p 

and q are polynomials in 


Then, it is 


easily verified that 
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in TR 
in 1R , 

and therefore v defined in (3.8) is the solution of (3.7). 
the proof. H 


(3.10) 


Xp - 


p = R 
xx 


Xq - 


q - S 
xx 


This completes 


To get the characteristic polynomial of the second derivative operator we need 
to substitute (3.8) into (3.6) and make use of (3.10) to get 

! atR(l) + a |£ ( 1 ,y ) ] + b[S(l) + a §£ (l,y)] = 0 

a[R(-l) - a (-1 ,w)] + b [ S ( — 1 ) - a (-l,y)] = 0. 

From now on we suppose that N is even (for N odd similar arguments can be 
applied), so that we have R(l) = S(l) = R(-l) = -S(-l) = 1 and 
p(x,y) = p(-x,y), q(x,y) = -q(-x,y). Hence, we can state: 

Theorem 3.1: The complex number X ^ 0 is an eigenvalue of (3.6) iff 

y = Y- satisfies 

(3.12) 2[1 + ct || (1 ,y ) ]• [1 + a || (l,y)] = 0. 

Proof : The left hand side of (3.12) is the determinant of (3.11). Since 

we are looking for nontrivial solution of (3.6), such determinant must vanish. H 
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Now, define 

g(u> - 1 + a (1 »y ) 

h(y) = 1 + a (l,y). 

It is not difficult to check that g and h are polynomials in y of 
N 

degree y . In order to show that the roots of g(y) and h(y) are 
real negative and distinct, we use the notion of a positive pair (see [5] and 
[8]). Two polynomials form a positive pair if their roots are real negative 
and interlaced. We shall prove, for instance, that g(y) and p(l,y)/y 
form a positive pair. To show that, we first need the following result (we 
recall that has been defined in Theorem 2.1). 

Lemma 3.2: Let 

2 

(3.14) f (y ) = g(y 2 ) + ay ■ ) ], 

y 

where g is defined in (3.13) and p in (3.9). Then f is a Hurwitz 
polynomial (i.e., all its roots lie in the left side of the imaginary axis) 
provided a is sufficiently large. 

Proof : By the definitions (3.9) and (3.13), one easily gets 

,, c / \ , . r D (2k+l) n x 2k+2 r c 2(k) / ..v 2k _ 

(3.15) f (y ) = 1 + a \ R (1) y + ay l R (1) y 

k=0 k=0 

-1+a I R (m) (l)y mfl . 
m=0 
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We show that f is the characteristic polynomial relative to the pseudo- 
spectral approximation of a hyperbolic problem. In fact, define 

p(x,y) = l R^ m \x)y IItfl , 
m=0 

then it is readily verified that 

(3.16) — p(x,y) = p (x,y) - ctp (1 ,y )R(x) , 

U x 

and that the roots of f (y ) * 1 + ctp ( 1 ,y ) = 0 give the corresponding eigen- 
values. Now, (3.16) actually is the eigenvalue problem associated with the 
hyperbolic equation 

3 Wj, 3 Wj, 

(3.17) a w N (l,t)R(x). 

With a proof similar to that of Theorem 2.1, where w^ plays the role of 
with g = 0 and Q = 0, it is possible to show that, for some norm II • II , 

j 2 

we have -2— llwjl <0 if a is suitably large. This implies that f is 
dt N 

Hurwitz. ■ 

As an immediate result of Lemma 3.2, we have the next theorem. 

Theorem 3.2: _If_ a is sufficiently large then the roots y of the 

polynomial g defined in (3.13) are real negative and distinct. 
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Proof : The theorem is a consequence of f being a Hurwitz polynomial. In 

fact, this is a necessary and sufficent conditon for g(p) and p(l,p)/p 
to form a positive pair (see [5], p. 228). In particular, the roots of g 
are real and negative. ■ 


In the same way, we can also prove: 


Theorem 3.3: If a is sufficiently large then the roots p of the 

polynomial h defined in (3,13) are real negative and distinct. 


Proof: It can be verified that the polynomials h(p) 

2 

a positive pair by showing that h(p^) + otp [■ 9 . ^ . ^* — -] 

P 

nomial. HI 


and q(l,p)/p form 
is a Hurwitz poly- 


Finally, by Theorems 3.1, 3.2, and 3.3, we can conclude with the following 
result. 


Theorem 3.4: _If_ a is sufficiently large then the eigenvalues X * 0 of 
the second derivative Chebyshev matrix with the boundary conditions (3.5), are 
real and negative. ■ 

It is easily verified that a turns out to be proportional to as Is 
also pointed out in [4], 
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4. ANALYSIS OP THE EIGENVALUES AND NUMERICAL EXPERIMENTS 

In this section, we analyze the behavior of the eigenvalues of the 
(N+l)x(N+l) matrix associated with the scheme (1.3). Applying the same proof 
of Theorem 2.1 in Section 2 to the equation (2.5) with g = 0, we get 



Figure 4.1 - Eigenvalues in the complex plane using scheme (1.3). 
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This says that such eigenvalues have negative real parts. In Figure 4.1, they 

o 

are plotted for N = 9 and a = N . The distribution in the complex plane 
is similar to that of the eigenvalues corresponding to the NxN matrix 

associated with the system (1.2). The extra eigenvalue coming from (1.3) is 
real, negative and its magnitude is proportional to N^. If R N (A) is the 
N 1 ** 1 degree characteristic polynomial related to (1.2) (see [2] for the 
explicit expression of the coefficients) it is easy to check that the eigen- 
values of (1.3) are the N+l roots of the equation 

(4.1) + aR N (X) = 0. 

The eigenfunction corresponding to the root \ of (4.1), up to a normaliz- 
ing constant, takes the following form 

(4.2) u(x) = l h^ k \x) X N k , 

k=0 

where h(x) = T'(x)(l+x). 

N 


To discretize in time (2.2), we can use the second order Runge-Kutta method. 
The analysis of the stability of the method, based on the knowledge of the 
eigenvalues of (1.3), gives an upper bound on the time step At. By choos- 
ing a proportional to N 2 , the restriction on At is, with a good 
approximation, represented by the formula 


At < 


• 4N 


(4.3) 


a - 
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2 2 
Therefore, by taking a = N , condition (4.3) says that At £ 3.3/N . 

Such restriction is slightly more severe than that obtained by exactly impos- 
ing the boundary condition in x = 1. Actually, in this last case, we had 
2 

At < 17/N (see [2]). The more restrictive condition on At is due to 

the presence of the real eigenvalue with the largest magnitude. One could 

think that this result negatively influences time discretization for scheme 
(2.2). Nevertheless, we argue that this is not the case. In fact, consider 
for instance problem (2.1), when the initial guess is f(x) = l-cos(x-l) 
and g = 0. We discretize the equation by collocation at the Chebyshev 
nodes x j , j = 1,***,N. Two different conditions are tested in x = 1, i.e. 


a ) Ujg( 1 , t ) — 0, 

(4.4) 

9u 3u 

b) g-jj— (l,t) = (l,t) - a u N (l,t). 


2 

We take N = 8, a = N and t£[0,T] with T=l, and we evaluate the error 
E as in (1.4) using both the schemes, respectively obtained by imposing con- 
ditions a or b in (4.4). Second order Runge-Kutta is used for time dis- 
cretization. Figure 4.2 shows the behavior of the error versus At. As the 
analysis of the eigenvalues pointed out, by increasing At using condi- 
tion b, instability occurs earlier than using condition a. 
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| Figure 4.2 - Comparison of the errors versus At using different boundary 

conditions. 

Nevertheless, the error relative to condition b is definitely lower than the 
other. Furthermore, the choice of a higher At is not appropriate because 
the time discretization error dominates. 

For the same example, Table 4.1 shows the error when T = 1 for various 
f choices of N. 

i 
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At = . 

01 

> 

ri- 

ll 

• 

001 

At = . 


N 

Condition a 

Condition b 

Condition a 

Condition b 

Condition a 

Condition b 

8 

. 1649E-02 

. 1 124E-02 

. 1644E-02 

. 1 122E-02 

. 1644E-02 

. 1 122E-02 

16 

. 2076E-03 

. 2039E-03 

. 1996E-03 

. 1962E-03 

. 1995E-03 

. 1961E-03 

32 

.6837E-04 

OVERFLOW 

.3846E-04 

.3653E-04 

.3842E-04 

.3649E-04 


Table 4.1 - Comparison of the errors for different At and N. 


Similarly, Table 4.2 shows the error when different values of T are used 
and At = .001. In almost all the cases the use of condition b is pre- 
ferred, especially when large values of T are considered. 



T = 

.5 

T = 

2 

T = 

"10 

N 

Condition a 

Condition b 

Condition a 

Condition b 

Condition a 

Condition b 

8 

. 1107E-02 

. 1200E-02 

.8904E-03 

.7875E-03 

. 2010E-07 

.4098E-09 

16 

.2594E-0e 

.2540E-03 

.9567E-04 

.8506E-04 

.8312E-10 

. 2718E-13 

32 

.4168E-04 

.3904E-04 

. 1 187E-04 

. 1174E-04 

.2231E-14 

.9872E-19 


Table 4.2 - Comparison of the errors for different T and N. 


Similar results can be obtained when time dependent boundary conditions are 
considered. 

We conclude this section by discussing preconditioning for the matrix 
corresponding to (1.3). For the matrix resulting from scheme (1.2) an 
efficient preconditioner was proposed in [3]. Such preconditioners can be 
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written as a product of two NxN matrices Z and D, where D is the up- 
wind finite-differences matrix at the collocation nodes and Z is a shift in 
the space of polynomials of degree N-l, from the values at the staggered grid 
points to the values at the initial grid. The eigenvalues after precondition- 
ing are real positive and between 1 and ^ • An analogous result holds 
for the (N+l)x(N+l) matrix corresponding to the scheme (1.3). As precondi- 

/» A A A* 

tioner for such matrix we set ZD, where Z and D take respectively 

the form 



The preconditioned eigenvalues can be explicitly computed also in this case. 
They are 


m sin 


X 0 1; A m 


2N 


. imr 

sia If 


, m = 1 , • • • ,N. 


In particular 1 < X < 5- . The corresponding eigenfunctions, up to a 

— m z 

multiplicative constant, are 


U m (x) ' V x) " 1 + a(X +0 ’ " ' 

m 


The preconditioner presented above is particularly suggested when steady state 
solutions of problem (2.2) have to be computed. 
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